Theoretical description of high-order harmonic generation in 

solids 



A. F. Kemper, 1 ' 2 'F] B. Moritz, 1,3 J. K. Freericks, 4 and T. P. Devereaux 1 ' 2 

1 Stanford Institute for Materials and Energy Science, 
SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA 
2 Geballe Laboratory for Advanced Materials, 
Stanford University, Stanford, CA 94305, USA 
3 Department of Physics and Astrophysics, 
University of North Dakota, Grand Forks, ND 58202, USA 
^Department of Physics, Georgetown University, Washington, DC 20057, USA 

(Dated: January 17, 2013) 

Abstract 

We consider several aspects of high-order harmonic generation in solids: the effects of elastic 
and inelastic scattering; varying pulse characteristics; and inclusion of material-specific parameters 
through a realistic band structure. We reproduce many observed characteristics of high har- 
monic generation experiments in solids including the formation of only odd harmonics in inversion- 
symmetric materials, and the nonlinear formation of high harmonics with increasing field. We find 
that the harmonic spectra are fairly robust against elastic and inelastic scattering. Furthermore, 
we find that the pulse characteristics play an important role in determining the harmonic spectra. 

PACS numbers: 
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I. INTRODUCTION 



The rapid development of time-resolved pump-probe experiments such as angle-resolved 
photoemission spectroscopy (ARPES) 13 , optical reflectivity 4 and X-ray scattering^ has 
expanded greatly the use of strong ultrashort laser pulses in the study of condensed mat- 
ter systems. The advent of these new experimental techniques requires theory to properly 
describe the non-equilibrium physics observed in these measurements. The strong tempo- 
ral dependence of these problems makes analysis in the frequency domain difficult which 
is more common for problems in equilibrium or within linear-response. A variety of meth- 
ods have been developed or adapted to describe the behavior of these systems directly in 
the time domain including semiclassical approximations,^^ Floquet theory^ and direct nu- 
merical integration of the time-dependent Schrodinger equation. 12 Each of these methods 
has limitations varying from lack of numerical feasibility to the inability to explicitly treat 
time-translation invariance breaking due to the pump and/or probe fields. 

The non-equilibrium Keldysh formalism, developed by Kadanoff and BaympI and 
Keldyshp31 can treat broken time-translation invariance explicitly making use of the so-called 
Keldysh contour (see Fig. [TJ. Quantities such as the single-particle Green's function are de- 
fined along this contour with multiple arguments that take appropriate real and imaginary 
times along the contour. Characteristics of the pump pulse, such as the driving frequency 
or pulse length and shape, as well as material specific tight-binding or more realistic band 
structure parameters for the equilibrium system can be included directly allowing for a 
closer match between theoretical calculations and realistic experimental conditions. The ab 
initio inclusion of nonlinear effects of the driving field, rather than treating it as a pertur- 
bation, is necessary for a proper description of the electronic system as it is driven out of 
equilibrium. This Keldysh formalism has been adapted and used to study the effect of elec- 
tronic correlations on systems out of equilibrium within dynamical mean-field theoryj^HSS 
continous-time Quantum Monte Carlo^I^ and cluster perturbation theory^!. A number of 
other methods exist to treat correlated and uncorrelated time-domain phenomena, includ- 
ing non-equilibrium exact diagonalizatiorP^ numerical renormalization group 26 and density 
matrix renormalization group with various benefits and limitations. 

An ideal application for this Keldysh approach is in the area of high-order harmonic 
generation (HHG) where one need only consider the non-linear effects of the pump, and 
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can neglect details associated with probe pulses and the often complicated evaluation of 
a proper pump-probe measurement cross-sections. High-order harmonics result from non- 
linear light-matter interactions and in the present example from the manner in which the 
sample current varies in time due to the pulse characteristics of the driving field. Used in 
atomic gases for the production of attosecond pulses,^ HHG opens avenues for time-domain 
experiments with high temporal resolution, and in particular are short enough to resolve 
electronic interactions which occur on these fundamental timescales. A second use for high- 
order harmonics is in the production of a wide spectrum of frequencies not available in the 
fundamental frequencies of standard lasers. By applying strong driving fields, the harmonics 
generated can be harvested for further use as a pump or probe. 

Trains of attosecond pulses that can then be used for further experiment ^ 30 1 31 1 can be 
generated through the interaction of a femtosecond laser pulse with a jet of gas. In this 
process the pulses of light are created in three-steps: (1) field-assisted tunnel ionization 
of an atom, (2) acceleration of the ionized electron back to the atom by the change in 
sign of the field as the driving pulse oscillates and (3) creation of a train of high-energy 
photons as the electron re-collides with the charged ion. 1 32 1 33 1 In addition to this generation 
in gaseous media, HHG can be accomplished through non-linear response in solids. The 
availability of ultra-short laser pulses with relatively strong fields has facilitated studies in 
solid state materials while mitigating the risk of sample damage. This capability has led 
to recent observations of HHG in ZnO crystals. 34 Their bulk crystalline nature makes HHG 
in these materials fundamentally different from what occurs in atoms or small molecules. 
The electron ionization in an atom is replaced by electron-hole pair excitation between the 
conduction and valence bands, respectively. The charged particles are then accelerated in 
the field with a non-linear response due to the Bloch oscillations which occur when the field 
strength is large enough, as discussed in some detail by several authorsP^ 

This is fundamentally different from the single-atom picture, where the harmonics are 
generated by the stimulated ionization and recombination of an electron in the outer shell of 
an atom. In a solid, this would correspond to electron-hole pair excitation and de-excitation. 
However, as we will show, high-order harmonics can be generated solely via the acceleration 
of the conduction electrons by the field, and the Bragg scattering that occurs repeatedly 
in high field.^ The end result is HHG of the fundamental frequency of the applied pump 
pulse. Although this has been previously studied for non-interacting systems J 9 * 10 * 34 ! using 
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the non-equilibrium Keldysh formalism, we can go beyond the semi-classical approach and 
consider the effects of elastic and inelastic scattering on the HHG spectra. 

A better understanding of HHG in solids requires a proper description of the underlying 
electronic states in the material. In atomic gases, these states have been modeled in a 
number of ways, from partial-wave expansion of the electron cloud to other descriptions 
utilizing first-principles methods.^! In solids, Bloch states are a basis in which one can 
describe the itinerant conduction band electrons and valence band holes. In the presence of 
a strong driving electric field, the states are modified Bloch states that, as the field energy 
becomes comparable to the bandwidth, become localized.^. The behavior covering the full 
range of fields, is captured by the inclusion of the field via the Peierls' substitution. 

The paper proceeds as follows. In Section II, we present relevant details of the non- 
equilibrium Keldysh formalism as it pertains to HHG. Section III discusses the effects of 
impurity and phonon scattering. In section IV we consider the effect of different pulse 
characteristics on HHG and compare and contrast the resulting HHG spectra for two exper- 
imentally relevant cases. Section V focuses on the time structure of the individual harmonics. 
Finally, section VI studies the HHG spectra from a more realistic, complex band structure 
appropriate for Si. 

II. METHOD 

We shall consider a system of electrons in the time-domain using the non-equilibrium 
Keldysh formalism. The appropriate expressions were derived in detaiP^ESl previously and 
we will present the relevant equations here for the specific case of a spatially uniform field. 
The system starts in an initial equilibrium state at t = t min , at an initial temperature T, 
and evolves until a final time t = t max . We describe the system on the Keldysh contour 
shown in Fig. [I] The time arguments of the double-time single-particle Green's functions 
and related quantities lie on the contour. Properly selecting the time arguments yields 
the time-ordered, anti-time-ordered, lesser and greater Green's functions, and by extension 
the retarded, advanced, and Keldysh Green's functions (see Ref. IT%]) . For example, the 
lesser Green's function (G < ) is defined with arguments t and if lying on the upper and 
lower branches of the contour, respectively. We use a particular implementation of the non- 
equilibrium Keldysh formalism where we discretize the Keldysh contour. The discretization 



4 



t 



mm 



t 



max: 



'mm 



-0 



Figure 1: The Keldysh contour used in the calculation of HHG spectra. 



is done in steps of 5 t and i5 T on the real and imaginary parts of the contour, respectively. This 
is different from the previous implementation by Turkowski and Freericks, who evaluated 
the temporal integrals for the Dyson equation^ directly.^ However, the direct evaluation of 
integrals is limited to the case of simple band structures, or, the infinite-dimensional case 
considered in the dynamical mean-field theory treatment of the Falicov-Kimball modelP^ 
Here, we numerically evaluate the integrals on a fine grid. 

The bare non-equilibrium Green's function Go under the influence of a driving field can 
be obtained through the Peierls' substitution^ in the standard definition of the double-time 
Green's function 



G {k;t,l/)=i[f j ;-e c {t,lf)] exp 



where f% is the Fermi function f% = ( e e ("/ T + 1 ) , T is the temperature (we work in 



i J <lt f I />• 

-l 



A(t) 



(1) 



units where the Boltzmann constant ks — 1), t and t' lie on the Keldysh contour, 9 is 
the contour-ordered version of the standard Heaviside function, e(k) is the single-particle 
dispersion, A(t) is the vector potential, which is related to the field by A(t) = — J E(t)dt, 
where we use e = c = h = 1 and work in the Hamiltonian gauge. The integral in the 
exponent runs along the Keldysh contour. 

To illustrate the influence of the electronic self-energy on HHG in solids, we include 
the effect of impurities and electron-phonon scattering. For impurities, we consider dilute 
point-like scatterers, for which the self-energy is®l 



£(M')=n^ 2 ^G (p; M '), 



(2) 
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where n« and V{ denote the impurity concentration and potential, respectively. For electron- 
phonon scattering in the Migdal limit of the Holstein model 41 , we consider a single optical 
phonon mode with frequency Vt coupled with a constant g, for which the self-energy is^ 

E(t, = ig 2 D (t, t') Go(p; t, (3) 



D (t,f) = -i 



[n B (0) + 1 - 9 c (t,t')]e iQ W + [ nB (0) + c (t,OK 



in(t-t') 



(4) 



where = l/(exp(f2/T) — 1). The full Green's function is then self-consistently calcu- 

lated through Dyson's equation expressed on the Keldysh contour as 

G{k; t, t') = G (k; t, t') + j) dt l( it 2 G {k; t, ti)E(t u t 2 )G{k; t 2 , t') (5) 

After discretization, this equation can be cast in the form of a matrix equation 

G = G + Go-(A tl SAi 2 )-G (6) 

where all times including those of the differentials A tl 2 lie on a branch of the Keldysh 
contour. This matrix equation can be solved using standard linear algebra techniques. 

The process by which high-order harmonics are generated in undoped semiconducting 
crystals begins with tunnel ionization of electrons from the valence band into the conduc- 
tion band by the pump laser. Next, this same pulse accelerates the electrons, creating an 
oscillating current. Here, the details of the tunnel ionization process are neglected and 
instead we focus on the motion of electrons only after they have been promoted to the con- 
duction band. In principle, there is a time dependence due to the tunneling process since the 
tunneling rate depends on the field, which is time dependent. However, since the tunneling 
rate depends exponentially on the field, the majority of electrons injected into the conduc- 
tion band will tunnel there when the field is largest, and will be essentially instantaneous.^ 
The subsequent driving force, in the absence of scattering will give rise to a similar harmonic 
spectrum as if the electrons had been there to begin with, albeit with somewhat broader 
harmonics and some higher frequency content due to the nearly instantaneous injection of 
the electrons into the valence band. Finally, as pointed out by Golde et al.j^the interaction 
between the excited electrons and the valence band holes leads to a contribution to the 
harmonic spectrum, in particular at very high harmonic orders. Here, we shall focus on the 
intra-band contribution, which is several orders of magnitude stronger. 



The instantaneous current j(t) created by the applied field is related to the lesser Green's 
function by 

j(t) = v(k-A(t))^mG<(k;t,t') (7) 

k 

— * — * — * 

where v(k) = de(k)/dk is the group velocity of a wave packet. We assume that the 
accelerated electrons radiate proportional to the Fourier transform of the instantaneous 
current, with the HHG spectrum proportional to 
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(8) 



We take the Fourier transform of the instantaneous current after interpolating with a shape- 
preserving Akima spline^. The self-energy and current are calculated on a 100 x 100 grid 
of momentum points in the two-dimensional square Brillouin zone, and 24 x 24 x 24 grid in 
the three-dimensional Brillouin zone appropriate for the Si band structure, sufficient to give 
converged results for the current. 

Since we are considering a laser pulse, a DC contribution to the electric field is excluded, 
as required by Maxwell's equations; any DC component does not propagate. Note that this 
restriction applies to the specific case of a propagating pulse; a DC contribution would be 
allowed for statically applied fields. A corollary to the absence of a DC component of the 
field is that the current is also zero at long times. This fact, combined with the relatively 
simple form of the current for a non-interacting system allows us to limit the length of the 
Keldysh contour to just the length of the pulse. Our discretization on the real-time branches 
of the Keldysh contour varies depending on the section of the paper, but is always kept above 
or equal to 5 t = O.OlSt^ 1 , allowing up to the 30 th harmonic to be visible (according to the 
Nyquist theorem)^. We always use a discretization of 10~ 2 T -1 along the imaginary-time 
spur, and set the simulation temperature to 0.1 eV. All times are measured in units of t^ , 
and applied fields in units of th- This can be converted to real units using appropriate values 
for the material under study. 

For sample values of the hopping th = 1 eV and material lattice constant or characteristic 
length scale ao = 5 A, a field strength of E = th corresponds to approximately 3.2 MV/cm 
(32 mV/A); the temporal unit in this work, t^ 1 would correspond to approximately 4.14 fs. 
Note that these values are not uniformly applicable; th and a need to be adjusted for the 
material under study. 



III. SCATTERING BY IMPURITIES AND PHONONS 



Motivated by the transition metal oxides, we consider a square lattice with a model band 
structure (in the absence of an electric field) given by 

e(k) = —2th(cos(k x ) + cos(k y )) + 4t' h cos(k x ) cos(k y ) — \i 

where t' h = 0.3t h and fi = —0.255th- We treat the system as isolated — unconnected to a 
bath — so the chemical potential remains at the initial equilibrium value. The filling mainly 
affects the amplitude of the current, and does not affect qualitatively the results . The 
system is driven by a Gaussian-shaped pump pulse centered at t — with a duration of 
4.24£-\ and a single underlying fundamental frequency u a = (9/5)tc which corresponds to 
roughly 9 full cycles under the Gaussian envelope. Since it is a pulse with finite width, the 
Fourier transform will contain a range of frequencies centered around u a , as will be discussed 
below. 

A. Harmonic generation from non-interacting electrons 

We calculate the current as a function of time for an oscillating field with a maximum 
field amplitude E max , as shown in Fig. [2] At the beginning of the pulse, the applied field 
amplitude is small and the field at that amplitude does not act very long; thus, the generated 
current is proportional to the field (up until t ~ —2t^), with equivalent current response 
exhibited at the end of the pulse. At larger amplitudes near the center of the pulse, the 
electrons are driven by the field for a sufficient period of time to hit the Brillouin zone 
boundary and Bragg scatter to the next zone, where their velocity points in the opposite 
direction. 

This phenomenon, commonly referred to as Bloch oscillations in condensed matter 
physics, originates from the periodicity of Bloch states. In the simple case of a constant 
DC field, the Peierls' substitutional becomes k — )■ k + Et and the 2ir periodicity of the 
band structure means that the frequency of Bloch oscillations satisfies ojb = 2-n/E. For 
oscillating fields, one can define an instantaneous Bloch frequency corresponding to the in- 
stantaneous field. Thus far, Bloch oscillations have not been observed directly in metals 
due to the experimental difficulty in resolving such short timescales and current degradation 
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Figure 2: a) Field profile A(t) for a multi-cycle Gaussian pulse with a maximum field amplitude 
of E max = 20th- (b) Instantaneous current response to this applied field. After the initial linear 
response, the field strength exceeds the point where Bloch oscillations occur, resulting in nonlinear 
behavior (see text), (c) Semi-log plot of the HHG spectrum (see text), showing peaks at odd 
harmonic frequencies. 

through impurity scattering, although oscillations have been observed in ultrasmall Joseph- 
son junctions,^ semiconductor superlattices, | 45 | 4S | and optically trapped atomsP^ However, 
the HHG spectrum is a direct consequence of the Bloch oscillations and provides an indirect 
measurement of them in solid-state systems, as discussed by various authors ! 9 I 1Q I 34 I 

As the field strength grows, the effects from Bloch oscillations intensify, as seen near 
the center of the pump pulse which displays several oscillations of the current in a region 
where the field does not change direction. A time-varying current due to the applied field 
leads to radiation, we thus perform a Fourier transform of the current to obtain the HHG 
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spectrum (shown in Fig. [2^c))P^ On a semi-logarithmic scale, we clearly observe peaks near 
odd multiples of the central frequency of the applied field (to a ). The fact that odd multiples 
dominate can be obtained from simple symmetry arguments and agrees with experimental 
observations. For inversion-symmetric systems, the velocity is anti-symmetric under parity, 
which leads to a cancellation of even harmonics of the current.^ The peaks have a finite 
width resulting from the finite pulse length where as a continuous pulse would lead to sharp 
peaks at exactly odd harmonic frequencies. 




Harmonic order 

Figure 3: The HHG spectrum at odd harmonics of the driving frequency for various field amplitudes 
using the 2D square lattice tight-binding model. 

Figure [3] shows the HHG spectrum at integer multiples of the driving frequency for various 
field amplitudes. For the smallest fields, the majority of the electrons do not Bragg scatter 
from the Brillouin zone boundary and there are neither signatures of Bloch oscillations nor 
high-order harmonics. This can be seen in Figure [3] by comparing the relative heights of the 
1 st and 3 rd harmonics; for the smallest field shown (E max = £/J the ratio of the 1 st to the 3 rd 
harmonic is roughly 10 -6 ; at this point, the current is entirely within the perturbative limit. 
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This ratio increases as a function of field with the largest fields showing a ratio approaching 
one. For the largest fields, the 1 st is not necessarily the largest, which can be understood 
from the time structure of the individual harmonics discussed below. 

B. Elastic scattering 

To study simple effects from elastic scatterers on HHG, we consider a pump pulse whose 
field amplitude (A max =20th) produces a sufficient number of odd harmonics so that variations 
in the HHG spectrum as a function of impurity scattering strength can be tracked reasonably 
well. Figure [4] shows the HHG spectra for increasing impurity strength, where we define a 



co/co 




Figure 4: Left: Current response for various impurity strengths, c = UiV- 2 , offset for clarity. Right: 
Corresponding HHG spectra, absolute (bottom) and normalized to the 1 harmonic (bottom) 

measure of impurity strength as c = n^V^ . As c increases, the major effect is an overall 
decrease in the magnitude of the instantaneous current, as can be seen in Fig.m Inaddit ion, 
near the center of the pulse where the Bloch oscillations are most noticeable, one observes a 
change in the current profile for the largest values of c. The harmonic spectra are affected in 
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a few ways. Most noticeably, the "noise floor" where individual harmonics can no longer be 
resolved (and the spectra go flat) lifts as the impurity scattering increases. This is due to the 
impurity scattering obscuring the smallest components of the electron motion (corresponding 
to the highest harmonics), and thus raising the background noise. Additionally, any small 
features in the clean spectra due to harmonics of the band structure are rapidly obscured.^ 

C. Inelastic scattering 
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Figure 5: Left: Current response for various electron-phonon coupling strengths, offset for clar- 
ity. Right: Corresponding HHG spectra, absolute (bottom) and normalized to the 1 st harmonic 
(bottom) 

In addition to impurity scattering, inelastic scattering can play a role in determining 
the current, and thus the HHG spectra. In particular, in semiconductors scattering due 
lattice vibrations (i.e. phonons) can play an important role.^Sl Here, we consider the effect 
of a single phonon mode on the generated spectra (as in the Holstein model) Due to 
computational restrictions, the pulse strength is limited to A max =5/f:/ l . We use a phonon 
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frequency of Q = 0.2th, and define the electron-phonon coupling strength as A = 2g 2 /QW, 
where W is the bandwidth and g is the electron-phonon coupling constant. With these 
parameters, we find the current and HHG spectra shown in Fig. [5] Similarly to the results 
for elastic scattering, although the overall current is strongly renormalized, the scaled HHG 
spectra only differ qualitatively at the highest frequencies, where the individual harmonics 
can no longer be resolved. 

IV. EFFECTS OF PULSE CHARACTERISTICS 





Figure 6: The field profile for the multi-cycle pulse (left) and single-cycle pulse (right). Inset: 
Fourier transform of the vector potential (see text). 




2 4 6 8 10 12 14 2 4 6 8 10 12 14 

^applied ^center 

Figure 7: The HHG spectra for a multi-cycle (left) and single-cycle (right) pulse. The color scale 
is logarithmic in intensity. Please note that the horizontal axis is scaled by the central frequency 
of the pulse in both cases. 
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The HHG spectra can be affected by the particulars of the applied pulse. To illustrate 
this, we consider two types of pulses: the multi-cycle pulse discussed above, and a single- 
cycle pulse similar to those used in THz experiments.^ For both pulses, although there is 
no DC component to the electric field, there is still a choice of the carrier-envelope phase. 
The results here have an arbitrarily chosen phase; we present a detailed dependence of the 
results of said phase in Section V. 

The field profiles for the two pulses considered are shown in Fig. [6j The inset shows the 
Fourier transform of the vector potential obtained from the field profile. Due to the relatively 
large number of oscillations in the multi-cycle pulse, the Fourier transform is strongly peaked 
at the driving frequency u a ~ 9tc/5. On the other hand, the single-cycle pulse is temporally 
short, leading to a large number of Fourier components with the largest contribution at 
uo « 4t h . 

Figure [7] shows the HHG spectrum as a function of field amplitude and frequency. There 
is a clear distinction between the spectra associated with the multi-cycle and single-cycle 
pulses. As discussed above, the current induced by the multi-cycle pulse responds in odd 
harmonics of the fundamental frequency which is relatively well defined as seen in the inset 
of Fig. [6j However, for the single-cycle pulse, the response is more complex and we find 
two distinct regimes. When the field amplitude is small, the electronic system responds 
perturbatively as in the small amplitude regime for the multi-cycle pulse. As such, the dom- 
inant contribution to the HHG spectrum comes from that frequency which gives the largest 
contribution to the field, Uq. However, this linear response is limited to the fundamental, 
at a field amplitude where no harmonics can yet be observed. As the field grows larger, the 
electrons undergo Bragg reflection as discussed above, which leads to non-linear response in 
the current . Since the Fourier spectrum of the single-cycle pulse is very broad, the peaks in 
the Fourier-transformed current do not simply occur at odd harmonics of the fundamental 
applied frequency, but rather at frequencies that depend on the field strength and that are 
not linearly spaced. Additionally, the peaks observed are much broader than those of the 
multi-cycle pulse, which reflects the broad Fourier profile of the single-cycle pulse. Lastly, 
we note that the single-cycle pulse has a clear asymmetry about the pulse center. We have 
repeated the calculations with a symmetrized pulse, and find no qualitative difference in the 
HHG spectra. 

In addition to the field frequency and amplitude, it is possible to experimentally control 
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Figure 8: HHG as a function of carrier-envelope phase for (a) the 9-cycle and (c) single-cycle pulse. 
The HHG spectra are multiplied by a finite offset for clarity. Panels (b) and (d) show the vector 
potential A(t) for the 9-cycle and single-cycle pulses, respectively. 

the carrier-envelope phase. This is directly addressable within the formalism presented, with 
the small caveat that there cannot be a DC component to the electric field. Thus, we choose 
to vary the phase in the vector potential: 

A(t) = A max sin(ut + 4>) exp ) ■ ( 9 ) 

Figure [8] shows the variation in the HHG spectra due to the carrier-envelope phase for both 
the 9-cycle and single-cycle pulse. The development of clear harmonic orders hinges on the 
applied field being sufficiently wide to clearly define the carrier frequency. Thus, the 9-cycle 
pulse does not show significant variation with the phase. However, the Fourier transform of 
the single-cycle pulse is very wide; the pulse is not long enough for a clear development of 
the carrier frequency. This causes a strong dependence on the carrier-envelope phase. 



V. TIME STRUCTURE OF THE INDIVIDUAL HARMONICS 

The HHG process in solids is due to the repeated Bragg scattering as the accelerated 
electrons hit the zone edge. This leads to a simple time structure in the generated harmonics 
(see also Refs. H0|34f48l) . Figure [9] shows the windowed Fourier transform (also known as 
a Gabor transform^!!) for the 9-cycle pulse. The low order harmonics are mainly generated 
in the low-field regions of the pulse. As the pulse grows stronger, the low order harmonics 
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Figure 9: a) Vector potential applied, b) Windowed Fourier transform (Gabor transform) of the 
current, showing the power spectrum of the n th harmonic as a function of time, c) The individual 
harmonic contribution to the current. 

are obscured by the higher orders, as the Bragg scattering occurs more rapidly. The higher 
orders are thus temporally localized to the center of the pulse. This can be clearly seen from 
both the time-dependent power spectra as well as the individual harmonics' contribution to 
the current. 

Through temporal localization one can understand the increase of harmonic amplitudes 
above the 1 st . The low order harmonics are generated near the beginning and end of the 
pulse, and their duration is inversely related to the maximum field strength. Thus, for large 
field strengths, the time during which the 1 st order harmonic has a strong contribution is 
smaller than the next orders. 
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VI. SILICON BAND STRUCTURE 
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Figure 10: HHG spectra at integer multiples of the carrier frequency in Si for fields along (100) 
(top) and (111) (bottom), for varying field strengths. The units for the field are all in eV. 



To investigate the effect of a realistic band structure on HHG, we have utilized a model for 
the conduction bands of Si as provided in some detail in Ref. [521 where we have placed the 
chemical potential just within the conduction bands (at the minimum near the X point). 



Fig. 10 shows the amplitude of the Fourier transformed current at odd multiples of u a 
(similar to the result shown in Fig. [3]). As in the case of the simple 2D band structure, the 
observed intensity increases with the field amplitude. Furthermore, since the bandwidth of 
the tight-binding model is different from the 2D model band structure, the field strength at 
which the nonlinearity onsets is similarly different. According to Ghimire et al.^, inclusion 
of more complex terms beyond first order in the band structure leads to an enhancement 
or reduction of the harmonic content in the HHG spectra. Here we consider the multi-cycle 
pulse discussed above, applied along the (111) and (100) directions of bulk Si. The band 
structure in both directions is sufficiently complex that, after averaging over the full Brillouin 
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zone, there is no singular qualitative difference which can be attributed to a particular feature 
of the band structure. There is an overall difference in the structure of the low field harmonic 
spectra, which appear to decay faster in the (100) direction compared to the (111) direction. 
The lack of a strong distinction in the spectra agrees with more detailed calculations done 



VII. SUMMARY 

In summary, we have considered several aspects of harmonic generation in solids. Har- 
monic generation from solids shows promise as a future tool as a tunable light source for 
further experiments. We find that the harmonic spectra are sensitive to the particular details 
of the driving pulse, which indicates a need to carefully control the parameters of said pulse. 
On the other hand, the harmonic spectra are fairly robust against the influence of scatter- 
ing, with the main feature being a decrease in the overall amplitude of the current. Finally, 
we find that for a realistic band structure, the HHG spectra do not depend sensitively on 
direction due to the complexity of the bands in any direction; however, this may not be the 
case in materials with simpler band structures. Overall, the most important message from 
this work is the very robust nature of the HHG response in solids, which does not seem to 
depend very strongly on the band structure, the orientation of the field, the scattering, or 
the pulse characteristics, except in the special cases detailed above. 

Acknowledgments 

AFK, BM and TPD were supported by the U.S. Department of Energy, Basic En- 
ergy Sciences, Materials Sciences and Engineering Division under contract No. DE-AC02- 
76SF00515. JKF was supported by the U.S. Department of Energy, Basic Energy Sciences, 
Materials Sciences and Engineering Division under contract No. DE-FG02-08ER46542 and 
by the McDevitt bequest at Georgetown University. The collaboration was supported by 
the U.S. Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering 
Division under contract Nos. DE-FG02-08ER46540 and de-sc0007091, This work was made 
possible by the resources of the National Energy Research Scientific Computing Center (via 
an Innovative and Novel Computational Impact on Theory and Experiment grant) which is 




18 



supported by the U.S. DOE, Office of Science, under Contract No. DE-AC02-05CH11231. 
We gratefully acknowledge discussions with P. S. Kirchmann, A. M. Lindenberg and D. A. 
Reis. 



* kemper@stanford.edu 

1 A. L. Cavalieri, N. Miiller, T. Uphues, V. S. Yakovlev, A. Baltuska, B. Horvath, B. Schmidt, 
L. Bliimel, R. Holzwarth, S. Hendel, et al., Nature 449, 1029 (2007). 

2 F. Schmitt, P. S. Kirchmann, U. Bovensiepen, R. G. Moore, L. Rettig, M. Krenz, J. H. Chu, 
N. Ru, L. Perfetti, D. H. Lu, et al., Science 321, 1649 (2008). 

3 L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, M. Wolf, H. Berger, S. Biermann, 
and A. Georges, New. J. Phys. 10, 053019 (2008). 

4 D. N. Basov, R. D. Averitt, D. van der Marel, M. Dressel, and K. Haule, Rev. Mod. Phys. 83, 
471 (2011). 

5 A. Cavalieri, C. Toth, C. Siders, J. Squier, F. Raksi, P. Forget, and J. Kieffer, Phys. Rev. Lett. 
87 (2001). 

6 A. Rousse, C. Rischel, and J.-C. Gauthier, Rev. Mod. Phys. 73, 17 (2001). 

7 T. Tritschler, O. D. Miicke, and M. Wegener, Phys. Rev. A 68, 033404 (2003). 

8 D. Golde, T. Meier, and S. W. Koch, Phys. Rev. B 77, 075330 (2008). 

9 D. Golde, M. Kira, T. Meier, and S. W. Koch, Physica Status Solidi (b) 248, 863 (2011). 

O. D. Miicke, Phys. Rev. B 84, 081202(R) (2011). 

1 F. H. M. Faisal and J. Z. Kamihski, Phys. Rev. A 56, 748 (1997). 

2 P. L. DeVries, J. Opt. Soc. Am. B 7, 517 (1990). 

3 L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics (W. A. Benjamin, Inc., NY, 1962). 

4 L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964). 

5 V. Turkowski and J. K. Freericks, Phys. Rev. B 71, 85104 (2005). 

6 J. K. Freericks, V. Turkowski, and V. Zlatic', Phys. Rev. Lett. 97, 266408 (2006). 

7 M. Eckstein and M. Kollar, Phys. Rev. Lett. 100, 120404 (2008). 

8 J. K. Freericks, Phys. Rev. B 77, 075109 (2008). 

9 M. Eckstein, A. Hackl, S. Kehrein, M. Kollar, M. Moeckel, P. Werner, and F. A. Wolf, Eur. 
Phys. J.: S.T. 180, 217 (2009). 

19 



20 B. Moritz, T. P. Devereaux, and J. K. Freericks, Phys. Rev. B 81, 165112 (2010). 

21 M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. B 81, 115131 (2010). 

22 B. Moritz, T. P. Devereaux, and J. K. Freericks, Comp. Phys. Comm. 182, 109 (2011). 

23 M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. Lett. 103, 056403 (2009). 

24 M. Balzer and M. Potthoff, Phys. Rev. B 83, 195132 (2011). 

25 T. J. Park and J. C. Light, J. Chem. Phys. 85, 5870 (1986). 

26 F. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 (2005). 

27 S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004). 

28 G. Alvarez, L. G. G. V. Dias da Silva, E. Ponce, and E. Dagotto, Phys. Rev. E 84, 056706 
(2011). 

29 P. Antoine, A. L'Huillier, and M. Lewenstein, Phys. Rev. Lett. 77, 1234 (1996). 

30 P. M. Paul, E. S. Toma, P. Breger, G. Mullot, F. Aug, P. Balcou, H. G. Muller, and P. Agostini, 
Science 292, 1689 (2001). 

31 M. Drescher, M. Hentschel, R. Kienberger, G. Tempea, C. Spielmann, G. A. Reider, P. B. 
Corkum, and F. Krausz, Science 291, 1923 (2001). 

32 J. L. Krause, K. J. Schafer, and K. C. Kulander, Phys. Rev. Lett. 68, 3535 (1992). 

33 M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L'huillier, and P. B. Corkum, Phys. Rev. A 49, 
2117 (1994). 

34 S. Ghimirc, A. D. Dichiara, E. Sistrunk, P. Agostini, L. F. Dimauro, and D. A. Reis, Nat. Phys. 
7, 138 (2011). 

35 G. Nenciu, Rev. Mod. Phys 63, 91 (1991). 

36 J. H. Davies and J. W. Wilkins, Phys. Rev. B 38, 1667 (1988). 

37 G. Mahan, Many-Particle Physics (Plenum Press, NY, 1990). 

38 L. M. Falicov and J. C. Kimball, Phys. Rev. Lett. 22, 997 (1969). 

39 R. Peierls, Zeitschrift fur Physik 80, 763 (1933). 

40 A. P. Jauho and J. W. Wilkins, Phys. Rev. B 29, 1919 (1984). 

41 T. Holstein, Annals of Physics 8, 325 (1959). 

42 H. Akima, J. ACM 17, 589 (1970). 

43 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, 
Numerical Recipes 3rd Edition: The Art of Scientific Computing (Cambridge University 
Press, New York, NY, USA, 2007), 3rd ed., ISBN 0521880688, 9780521880688. 

20 



L. Kuzmin and D. Haviland, Phys. Rev. Lett. 67, 2890 (1991). 

C. Waschke, H. Roskos, R. Schwedler, K. Leo, H. Kurz, and K. Kohler, Phys. Rev. Lett. 70, 
3319 (1993). 

V. Lyssenko, G. Valusis, F. Loser, and T. Hasche, Phys. Rev. Lett. 79, 301 (1997). 

M. Ben Dahan, E. Peik, J. Reichel, Y. Castin, and C. Salomon, Phys. Rev. Lett. 76, 4508 

(1996). 

S. Ghimire, A. D. DiChiara, E. Sistrunk, G. Ndabashimiye, U. B. Szafruga, A. Mohammad, 
P. Agostini, L. F. DiMauro, and D. A. Reis, Phys. Rev. A 85, 043836 (2012). 
V. Zhukov and E. Chulkov, J. Phys. Cond. Matter 22, 435802 (2010). 
A. Lindenberg, private communication (2011). 

K. Grochenig, Foundations of Time-Frequency Analysis (Birkhauser, 2001). 
P. Vogl, Journal of Physics and Chemistry of Solids 44, 365 (1983). 

J. K. Freericks, A. Y. Liu, A. F. Kemper, and T. P. Devereaux, Physica Scripta Volume T 151, 
014062 (2012), 1212.1610. 



21 



